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Abstract 


We propose a combined reconstruction-classification method for simultaneously recovering absorption and scattering in turbid 
media from images of absorbed optical energy. This method exploits knowledge that optical parameters are determined by a 
limited number of classes to iteratively improve their estimate. Numerical experiments show that the proposed approach allows for 
accurate recovery of absorption and scattering in 2 and 3 dimensions, and delivers superior image quality with respect to traditional 
reconstruction-only approaches. 
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1. Introduction 

Photoacoustic tomography (PAT) is an emerging technique 
for in vivo imaging of soft biological tissue [1]. This hybrid 
modality uses ultrasound to detect optical contrast, combin¬ 
ing the high resolution of acoustic methods with the spectro¬ 
scopic capability of optical imaging. To generate a PA im¬ 
age, a short laser pulse is shone into the object, the ultrasonic 
waves emitted following the heating of the tissue are measured, 
and an image of the absorbed optical energy field is recovered. 
Whereas purely optical methods suffer from poor spatial resolu¬ 
tion, acoustic waves propagate with minimal scattering and PAT 
can achieve 100 micron resolution at depths of several centime¬ 
tres. However, PA images provide only qualitative information 
about the tissue, and are not directly related to tissue morphol¬ 
ogy and functionality. The principal difficulty is that the PA 
image is the product of both the optical absorption coefficient 
(which is directly related to underlying tissue composition) and 
the light distribution (which is not). This severely restricts the 
range of applications for which PAT is suitable. 

Quantitative photoacoustic tomography (QPAT) aims to 
provide clinically valuable images of the optical absorption 
and scattering coefficients, or chromophore (light-absorbing 
molecules) concentrations from conventional PA images via an 
image reconstruction method [2]. A model of light propagation 
is required to relate the absorbed optical energy to the light fiu- 
ence and tissue parameters. The primary challenge of QPAT is 
solving the non-linear imaging problem. In particular, recover¬ 
ing the scattering coefficient is especially difficult to due to its 
weak dependence on the absorbed energy density. 

In this paper, we develop a method for solving the image 
reconstruction problem for QPAT by alternating reconstruction 
and segmentation steps in an automated iterative process. We 
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introduce a probabilistic model that describes optical properties 
in terms of a limited number of optically distinct classes, which 
may correspond to tissues or chromophores. These are identi¬ 
fied and characterized by a classification, or segmentation, algo¬ 
rithm. This approach allows for the use of information retrieved 
by the classification in the reconstruction stage, and vice versa. 
The aim of the reconstruction is to choose solutions for which 
the image parameters take values close to a finite set of discrete 
points. The aim of the classification algorithm is to progres¬ 
sively improve the parametric optical model, and correct for er¬ 
rors in the initial assumptions. Multinomial models have been 
employed previously in the related fields Diffuse Optical To¬ 
mography [3] and Electrical Impedance Tomography [4]. For 
QPAT, the main advantage is that this approach enables accu¬ 
rate recovery of both the absorption and scattering coefficients, 
simultaneously. 

2. Numerical methods 

2.7. Quantitative photoacoustic imaging 

A conventional PAT image is proportional to the absorbed 
optical energy 

H{r) = t{r)pa{r)fy(pa (r)) r g Q, (1) 

where r is a position vector within the domain Q, pa and p'^ 
are the optical absorption and reduced scattering coefficients, 0 
is the optical fiuence, and T is the Griineisen parameter. The 
Griineisen parameter represents the efficiency with which the 
tissue converts heat into acoustic pressure, and is often taken 
to be constant T(r) = 1, Vr g Q. The fiuence is dependent on 
the optical parameters and illumination pattern in the whole do¬ 
main. The problem of recovering the optical parameters (pa,p's) 
from a conventional PAT image is known as the quantitative 
problem. The optical absorption pa is of particular interest be¬ 
cause it is fundamentally related to underlying tissue physiol¬ 
ogy and functionality, and encodes clinically useful information 
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such as tissue oxygenation levels and chromophore concentra¬ 
tions. Conversely, the absorbed energy density H depends non- 
trivially on optical absorption, thus is not directly related to tis¬ 
sue morphology because it is distorted, structurally and spec¬ 
trally, by the non-uniform light fluence. 

2.2. The diffusion model of light transport 

In order to recover the optical parameters a model 

of light propagation within the tissue is required. For highly 
scattering media and far from boundaries and sources, a low or¬ 
der spherical harmonic approximation to the radiative transfer 
equation is suitable. The diffusion approximation is given by 

[5] 

(/i,-V.^(r)V)0(r) = ^(r), (2) 

where q{r) is an isotropic source term, and k = 1/3//' is the 
diffusion coefficient. 

We set Robin boundary conditions 

0(r) -r ^K{r)h • V0(r) = 0 r e dQ (3) 

where A accounts for the refractive index mismatch at the 
boundary. 


Substituting into the the objective function (4) leads to the dis¬ 
crete form of the objective function 

^=\ - {%, \ - {'Vj, I^a<py. 

i i 

(7) 

If a single illumination source is used and both absorption 
and scattering are undetermined, the problem is ill posed [2]. 
In this study, the non-uniqueness of the solution was removed 
by using multiple illumination patterns [6], thus the objective 
function must be summed over the number of sources. In the 
following, we have omitted this sum for ease of notation. Prior 
information regarding the solution can be included by adding a 
regularization term 

(8) 

j 

In the Bayesian framework, an image is obtained by maximiz¬ 
ing the posterior probability of the parameters, given the data: 

p(//a,//;|J'”) cx (9) 

Under this interpretation, the regularization term R is given by 
the negative log of the prior probability distribution 


2.3. Minimization-based QPAT imaging 


R{pa,p's) = -\ogp{pa,p[). ( 10 ) 


In this paper, we adopt a gradient-based minimization ap¬ 
proach to image reconstruction. Typically, both pa and //' are 
unknown and need to be recovered simultaneously from the ab¬ 
sorbed energy density. An objective function is defined, which 
measures the distance between the conventional PAT image 
and the data predicted by the model for the current estimates 

6=1 f (H'" - ( 4 ) 

^ Jq 

In order to treat the problem for a generic geometry, the Fi¬ 
nite Element Method (FEM) is employed, whereby a weak for¬ 
mulation of the diffusion approximation (2) is considered. A 
discretization of the domain is defined, and the fluence and op¬ 
tical parameters are expressed in terms of piecewise linear basis 
functions M,(r); ZiXiUi(r) forx e where are 

nodal coefficients and i = 1,..., A. 

We assume that the data is the absorbed energy density 
projected onto a particular basis j'Pyj, 

d'^ = [dJ,j=l,...,N), (5) 

dj = f H"’(r)'¥j(r)dQ = ('P, //“). (6) 

Jq 


2.4. Gradient calculations 

Cox et al. [7] have shown that, for the continuous case, the 
gradient of (4) with respect to Pa at position is given by 

^ =-(P{H"'-H)\^ + (11) 


where 0* is the adjoint light field. In the following, we derive 
the expression for the gradient in the discrete case. 

The sampled forward model can be expressed as a vector 

H=[HjJ=\,...,N\ 

Hj= f H(r)'P/r)dn = ('P, //'"), 

Jq 

= y/^aA f ^j(r)ui(r)u,(r)dQ = ^^ajii„, (12) 
IT dn 


where is a sparse matrix with entries /, k where the support 
of the basis functions 'Py(r), ufr), Uk(r) overlap. Taking the 
derivative of (7) with respect to pa , , we have 



(13) 


Choices for j'Pyj include: 

1. Point sampling 'Fy(r) = 6{r - rj), 

2. Piecewise-linear sampling 'Ey = uj, 

3. Sine sampling 'Ey = sinc(|r - rj\). 


Using the expression for the absorbed energy density (12), 


dPai 


eJC^(l>^pJC^ 


d(f> 

dpai' 


(14) 


where et is a vector of zeros with a single 1 in position i. Sub- 
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stituting into (13) gives 


(15) 


Substituting into (15) gives the expression for the derivative 
with respect to ^ 


dS 


■ VieJ CV + - Hj). 


The first term in equation (15) is 


ejCWJ - Hj) = eiCi<Pu{d'J - Hj) 

j j,i,k 

= J]MdJ-Hj) f '¥j(r)ui(r)udr)dQ 

j,k 

= ifi^E‘(d"’ - H) (16) 

where E' is given by a reordering of 


E[. = ^^^j{r)ui{r)u^{r)d£l. (17) 

Note that while is symmetric, in general is not. 

It remains to determine The discrete form of the DA 

HIda i 

model ( 2 ) assumes the form [ 8 ]: 


(M + K + F)0 = e, (18) 


88 

dHai 




(28) 


The derivative with respect to yu' ■ can be derived analogously: 


dS 

dfi'si 




(29) 


where 

= UiVuj -VukdQ, (30) 

and ^ = -l/3yu'^. Note that calculation of the gradient only 
requires two runs of the forward model. The forward problem 
was solved using the Toast++ software package [ 8 ]. 

Choosing point-sampling ^j(r) = d(r - rj), gives simply 
= I. In this study, we chose piecewise-linear sam¬ 
pling = Uj, so we had = E^ = and 


88 

dfiai 




(31) 


where 


^jk ~ ^ i ^ UiUjUj^dFl, 

(19) 

^jk — ^ 1 UiWuj • Vi/^dD, 

i 

( 20 ) 


( 21 ) 

Qj=y\di 1 UiUjdSl. 

V Jn 

( 22 ) 


i 


Taking the derivative of equation (18) with respect to the ith 
coefficient of fia, 

(M + K + F)^ =-V;^^ (23) 

Ofia i 

where 

K.,jk= f UiUjUkd^i (24) 

is given by the derivative of the system matrix. We define the 
adjoint field 0 * as the solution to the equation 

(M + K + F)(fi* = Q* (25) 

where 

Q* = Y,t^JC\d'J-Hj) (26) 

j 

is the adjoint source. Taking 0*- (23) - ^ • (25) we obtain 

J^n/d^(dJ - Hj) = (27) 


3. Reconstruction-classification method for QPAT 

A reconstruction-classification scheme is devised, which en¬ 
ables the recovery /i^ and yu' by approaching the image recon¬ 
struction and segmentation problems simultaneously. At each 
reconstruction step, we minimize a regularized objective func¬ 
tion, where the regularization term is given by a mixture model. 
At each classification step, the result of the previous reconstruc¬ 
tion step is employed to update the class parameters for the 
multinomial model. We alternate between reconstruction and 
classification steps for a fixed number of iterations. 


3.1. Mixture model for yu^ and yu' 

In this section we introduce a probability model for jUa and 
yu', which encodes prior knowledge about the optical parame¬ 
ters and allows us to bias the solution of the imaging problem 
accordingly. We assume that an array of labels can be deter¬ 
mined for each node, such that 


= 


1 if the ith node is assigned to the 7 th class; 

0 otherwise. ^ 


The labels constitute hidden variables on which the image pa¬ 
rameters are dependant. For each class j = 1,...,/, a mean 
vector nij = {jia ^ is defined, and the closeness of the 

optical parameters to the mean values is described by a covari¬ 
ance matrix 11 j e 

We assume that if = 1, the probability distribution for 
Xi = (fta vl^'s /) is given by a multivariate Gaussian distribution 

= (33) 


where 6j indicates the set of class parameters {mj, Hj). 
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The prior probability distribution of the class properties Oj is 
given by the conjugate prior to the Gaussian distribution. Prior 
information about the distribution of the class means or covari¬ 
ances can be encoded by choosing the parameters of the con¬ 
jugate prior accordingly. Using a non-informative prior for the 
class means we have oc 1. The conjugate prior distribu¬ 

tion for the covariance of a normal distribution is given by the 
normal inverse Wishart distribution 


NIW(y^-,r^-) 


-(y+J+l)/2 


exp 


-^Tr(r,. 2 ;Ti) 


(34) 


where d is the dimension of the domain, vj indicates the number 
of degrees of freedom, and T^ is a scaling matrix. If the prior 
is non-informative, then Vj = 0 and T^ = 0, and the probability 
distribution of the class parameters becomes 


P(^i) “ pil 


-W+l)/2 


(35) 


which is known as Jeffreys prior. 

The probability that the set of labels = {f/i,..., ..., is 

assigned to the ith node is given by a multinomial distribution 

p((i\A)=Y]Af. (36) 

j 

where Aj is the overall probability that a node is assigned to the 
jth class. Therefore the joint probability for (jc/, is given by 
the product 


is obtained by the fixing the labels to the maximum a posteriori 
estimate, given the results of the previous iteration 

MAP(^) = argmaxp(^|y-‘,6»'-\/-‘), (42) 

which is calculated in the classification step (see section 3.1.2). 
The weighting matrix is the Cholesky decomposition of 
where e ]^ 2 vx 2 v ^ sparse matrix of which the ith 

2x2 block along the diagonal is Xf if the ith element belongs 
to the fth class. 

In order to sphere the solution space, that is to render the 
space dimensionless, we performed a change of variables Pa 
Pa/l^ao and//' ^ 

for the optical parameters (in this study, we initialized to the ho¬ 
mogeneous background). Given the size of problem, we chose a 
gradient-based optimization method in order to reduce memory 
use and computational expense [9] . The minimization was per¬ 
formed using the limited-memory Broyden-Fletcher-Goldfarb- 
Shanno (L-BFGS) method [10], with a storage memory of 6 
iterations. 

3.1.2. Classification 

The purpose of the classification step is to update the multi¬ 
nomial model, using the result of the previous reconstruction 
step. First, the expected values of the labels are computed 
for the current class parameters (0\ A^) and image = {jiJ, p'J) 
(E-step). Then the model parameters are updated by maximiz¬ 
ing the posterior probability (M-step) 


p(x„ Q0, A) = p(x,|^,., 0)P(^,|4) = Y] YjP(Xi\0j)f ■ (37) 

j 


P(6/, AW) oc p{W\e, A)p(6, A). (43) 


By marginalizing over all possible values of the indicator vari¬ 
ables ^ij, a mixture of Gaussians model for the optical parame¬ 
ters is obtained 

p(xt\0,A)= f p(x,-,^,|^,i)d^-yT,-p(x,|^;). (38) 

J 

Finally, for independent nodes the prior of the image is given 
by 

p(x|6»,/t) = ns AjP(Xil0j). (39) 

i j 


E-step: 

The responsibility r\. is a measure of the probability that the ith 
node is assigned to the jth class. Using Bayes’ theorem and the 
Gaussian mixture model (38) we have 




V(Xi\^ij = l,0‘)p(^ij = 1 ) 

P(a:,|0, A) 

Ay(x'i0‘) ^ 


ZjA‘p(x‘J0‘) 

The expectation for the indicator values is 


nj’ 


(44) 


3.1.1. Reconstruction step 

The objective function takes the form of equation (8), where 
at iteration t of the reconstruction-classification algorithm the 
regularization is given by (equations (10) and (39)) 


E(i^ij\x\,ff,X) = J ^ijpd^ij = l\x‘,&,A‘)d^ij 

= 0 X pi^ij = 0\x\, 0, /) + 1 X pi^ij = l\x\, 0, A‘) 


Therefore the MAP estimate for the labels is 


(45) 


'R‘(Ma,n:) = -\ogp(x\0,A‘) = Y\\U(x - x)||2, 

where r is a regularization parameter and 


(40) 


^t+i 


1 if H. is maximum V j, 

U 

0 otherwise. 


which can be used in equation (42). 


(46) 


Xi = 


^ ^ij ■ 


my € 


(41) M-step: 

The parameters {6, X) are chosen in order to maximize the log 


MAP(^) 
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posterior 

= arg max log p(jc^|^, i) + log p(^, X). (47) 

{0,A) 

Averaging over all possible values of ^ gives 

log p(x‘\0, /t) + log p{0, ^) = J' log + log p{0, i) 

^ (48) 

Using Jensen's inequality [11] and ignoring terms which do not 
depend on {6, i), we obtain a lower bound for the log-prior 


3.3. Visualization of the results 

Results obtained using the reconstruction-classification 
method are displayed alongside scatter plots of the nodal val¬ 
ues recovered in the 2D feature space (//a,yu') (for example, 
see figure final column in 4). The positions of the class means 
nij = are identified by a cross, and the class covari¬ 

ances are represented by ellipses. These are colour coded by 
class, and are indicative of the clustering of image nodal values 
around the class means. 

4. Results 


^(6>, i) = sz Aj log (4yP(a-„|6»y)) + log p(i) + log p(61) 


I J 




log(/l;) + logCIEjl) - ^(jc,(„) - OTjO'E - mj) 


Z 


(a,-l)log(^,)- 


V; -h d 1 


logP,- 


(49) 


Maximizing !B(0, A) for Z; dj = 1 and using non-informative 
priors, we obtain the update rules for the model parameters 


4?' = 


yd. 


J N ’ 


W +1 

ntj = 


yd. ’ 

Zj( 'i 


^r+l 


Z; d.fXi - nijXxi - ttijf + r j 
I.iAj + vj + d + \ 


(50) 

(51) 

(52) 


3.2. Class means initialization 


The number of classes J and the class means ntj were ini¬ 
tialized by automatically segmenting the result of the first re¬ 
construction step and averaging over the segmented areas. To 
segment the image (for example, see figure la) we looked at a 
binned histogram of the image of and chose the value jUa ^ for 
which the number of occurrences was highest (figure Ic, col¬ 
umn 1). We found the first node index h for which the value jUa ^ 
occurs, and identified the corresponding scattering value yu' 
Having chosen a covariance matrix Yh, we computed a map of 
the multivariate normal probability of the images, with 

mean (qia i) (figure Ic, column 2). Then we selected a tol¬ 
erance level io\h at which to truncate the probability map, and 
selected all nodes with probability higher than the tolerance as 
belonging to the same class as node h (figure Ic, column 3). 
We repeated this process on the remaining nodes until all nodes 
were classified. Thus the number of classes was set to the num¬ 
ber of iterations, and the average of the optical parameters over 
each class was used to initialize the class means (figure lb). 


4.1. 2D validation and reconstruction 

We chose a numerical phantom defined on a 2D circular 
mesh with 1331 nodes and radius 25 mm. Four illumination 
sources were placed on the boundary at angles 0, 7r/2, n and 
37r/2 rad. In all cases the illumination profile was a normal¬ 
ized Gaussian with radius (distance from the centre at which 
the profile drops to Ije) 6 mm. The background optical pa¬ 
rameters were set to jUa = 0.01 mm“^ and yu' = 1 mm“^ Two 
circular perturbations of radius 6 mm were added in positions 
(6 mm, 10 mm) and (-6 mm,-10 mm) (figure 3a). The values 
of the perturbations were fia = 0.02mm“^ yu' = 1.5 mm“^ and 
jUa = 0.03mm“\yu' = 1.25 mm“\ respectively. The absorbed 
energy field was simulated for each illumination and 1 % white 
Gaussian noise was added (figure 3b). The class covariances 
were initialized to 


/ 10-6 0 
\ 0 10-1 


V7=1,...,3, 


(53) 


where the first variable was the absorption and the second was 
the reduced scattering. The parameters of the Jeffreys prior 
were set to T^ = Yj Vj, y(l) = 1 for the background class, 
and y(2, 3) = 10 for the perturbation classes. The number of 
classes and optical parameters were initialized using the class 
means initialization method (section 3.2) with tolh = 10-^ and 
Yh = Yj (53), and the labels were initialized to 1 for the back¬ 
ground class and zero for all other classes. The tolerance of the 
L-BFGS algorithm was set to tol = lO-^ and the total number 
of reconstruction-classification iterations was set to Maxit =10 
(figure 4). The regularization parameter r = 10-was cho¬ 
sen by inspection. For comparison, images were reconstructed 
without introducing a prior (figure 5); the images were recon¬ 
structed by minimizing (7) using the L-BFGS method with 
tol = 10-12. 


4.2. 3D validation and reconstruction 

We chose a 3D phantom analogous to the 2D case, de¬ 
fined on a cylinder with 27084 nodes, radius 25 mm and height 
25 mm. Two spherical inclusions of radius 6 mm were placed in 
(6 mm, 10 mm, 0 mm) and (-6 mm, -10 mm, 0 mm) (figure 6a). 
Illuminations sources were Gaussian in the xy-plane constant 
in the z-axis, with radius 6 mm and length 25 mm (figures 6b, 
6c). PAT images were simulated for 4 illuminations at the car¬ 
dinal points, and 1% noise was added to the absorbed energy 
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(a) 



Pa p's 

1 

0.009 

1.03 

2 

0.024 

0.15 

3 

0.017 

1.57 


(b) 




(c) 

Figure 1: Class initialization example: (a) original image of /da to which we apply the segmentation; (b) result of taking average image values over the segmented 
areas (c) first column, histogram of occurrences of values of jda in the portion of the image requiring segmentation - value with highest number of occurrences is 
/dah (indicated by a red cross); second column, probability density function with meaniju^/j,//'^) and covariance E/j; third column, labels identifying nodes with 
probability density higher than tolerance value tolh', each row corresponds to an iteration and a distinct class, so in this case J = 3. 


Set Maxit, tol 

Initialize optical and class parameters x, 6, A 
Initialize iteration count t = I and regularization term "RP = 0 
(no regularization) 

repeat 

Reconstruction: 

Update x\ minimize (8) using L-BFGS until 6 < tol 

if f = 1 then 

Initialize class means nij (section 3.2) 

end if 

Classification: 

E-step; compute expected labels ^ (44) 

M-step; update class parameters {6, X) (50)(51)(52) 
Update regularization term R! (40) 

Re-initialize jc = Jc (41) 
t <— f -f 1 
until t > Maxit 


Figure 2: Reconstruction-classification algorithm outline 


(figure 6d). The optical, covariance and reconstruction param¬ 
eters were set to the same values used in the 2D case. The class 
initialization parameters were set to tol/^ = 10“^ and 
Images were reconstructed by performing 10 iterations of the 
reconstruction-classification method (figure 7). 

5. Discussion 

5.7. Summary of findings 

We applied the proposed reconstruction-classification algo¬ 
rithm to a 2D numerical phantom with 3 tissues, a background 
and 2 perturbations (figure 3). The optical absorption was 
recovered reliably within a small number of iterations, and 
the scattering was recovered with sufficient accuracy after ap¬ 
proximately 10 iterations (figure 4). We compared the optical 
model with images obtained by the reconstruction-classification 
method, and by a traditional reconstruction-only (no regular¬ 
ization) method (figure 5). We found that the reconstruction- 
classification method delivered superior image quality, partic¬ 
ularly with regards to the scattering parameter. We applied 
the reconstruction-classification algorithm to a much larger 3D 
problem (figure 6) and observed similar results (figure 7) as in 
the 2D case. 

5.2. Choice of parameters 

The parametric optical model and classification algorithm in¬ 
troduce a number of parameters which require tuning by the 
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Figure 3: 2D model: (a) circular mesh, (b) absorbed energy for each illumination pattern. 
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Figure 4: 2D reconstruction-classification results at iteration 1 (first row), 5 (second row) and 10 (third row). Reconstructed values of jUa and //' (first and second 
column), labels recovered for perturbation classes (third and fourth columns), and scatter plot (fifth column). 


user. In addition to the regularization parameter, the parame¬ 
ters of the Jeffreys prior F and v and the initial guess of the 
class variances must be set before performing the classifi¬ 
cation. However, their significance is fairly intuitive, and with 
experience of a certain type of problem the choice of parame¬ 
ters becomes natural. Visualizing the class covariance matrix 
j as an ellipse, changing the value of F varies its eccentric¬ 
ity, and changing v varies the length of its axes. Further, given 
that in the first iteration the optical absorption is recovered with 
superior accuracy than the scattering, it is preferable to initial¬ 
ize the variance of the former to a smaller value than the latter, 
indicating greater confidence in the imaging solution. 

5.3. Initialization of the class means 

The purpose of the means initialization scheme is to increase 
automation of the method, so that minimum user intervention 


and no prior knowledge of the number of tissues or their op¬ 
tical properties is required. The algorithm simply performs a 
segmentation of the image, and then takes averages over the 
segmented areas to initialize the class properties (figure 1). Al¬ 
ternative segmentation techniques could have been employed, 
however the advantage of the proposed approach is that it di¬ 
rectly exploits the mixture of Gaussians model to identify the 
tissues. Our choice to investigate a node h with jda belonging 
to the bin with maximum number of occurrences leads to the 
background tissue being identified first, followed by the per¬ 
turbation tissues. The choice of the node index h could have 
been randomized, so that tissues are identified in random or¬ 
der. This approach is equally valid, however we found that in 
cases where tissue values were close together (such as after a 
single reconstruction-classification iteration) it was preferable 
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Figure 5: 2D model and reconstruction: first column, model of jds and ; second colunm: reconstructed values of jda and without multinomial prior; third 
column: reconstructed values of fia and with multinomial prior. 
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Figure 6: 3D model: (a) numerical phantom and perturbation locations, (b) all illumination sources, (c) cross section of optical parameters used to simulate the data 
for z = 0, (d) cross section of absorbed energy for each illumination pattern. 


to identify the largest classes first because the mean was esti¬ 
mated with greater accuracy for the classes with a larger num¬ 
ber of samples. Further, for a given image and tolerance level, 
our choice renders the result of the segmentation process unique 
and reproducible. 

5.4. Recovery of the scattering 

From the comparison with the reconstruction-only case with 
no regularization (figure 5), it is evident that the introduction 
of the parametric prior enables better recovery of the scatter¬ 
ing. The inconsistency between the quality of the recovered ab¬ 
sorption and scattering parameters in the non-regularized case 


is due to the weaker dependence of the latter on the absorbed 
energy density with respect to the former. This results in the 
scattering gradient being approximately an order of magnitude 
smaller than the absorption gradient. Although the problem can 
be mitigated by sphering the solution space, variations in the 
data due to the scattering often fall below the noise fioor. In 
the reconstruction-classification case, typically the absorption 
is recovered with good accuracy within a small number of iter¬ 
ations. Thus, the absorption takes values very close to the class 
means (resulting in small clusters), and the variance along the 
Ida direction converges to a small value. Given that the regular- 
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Figure 7: 3D reconstruction-classification results at iteration 1 (first row), 5 (second row) and 10 (third row). Reconstructed values of fia and //' (first and second 
column), labels recovered for perturbation classes (third column), and scatter plot (fourth column). 


ization term is weighted by the inverse of the covariance matrix, 
the dependence of the absorption gradient on the data becomes 
weaker at each iteration, until its magnitude is comparable or 
smaller to that of the scattering. In the iterations that follow, the 
descent of the data term of the objective function is primarily 
due to updating the scattering, which converges to the correct 
values. 

5.5. Computational demands 

Computational performance was found to be strongly depen¬ 
dent on the problem size. In the 2D case with 1331 nodes (fig¬ 
ure 4), the total reconstruction time (10 outer reconstruction- 
classification iterations) using Matlab on a 16 core PC with 
128GiB RAM was only 77 seconds. In the 3D case with 27084 
nodes (figure 7), the total reconstruction time increased linearly 
with the number of nodes, and on the same workstation was ap¬ 
proximately 3.7 hours. The increase in computation time was 
mostly due to much longer processing times for the L-BFGS 
algorithm in the reconstruction step. 

5.(5. Experimental application 

In experimental situations, prior information on tissue prop¬ 
erties may be held, such as knowledge of the characteristic opti¬ 
cal absorption and scattering spectra of chromophores of inter¬ 
est. These may be obtained from the literature [12], or gained 


through tissue sample measurements. This information could 
be used in one of two ways. Firstly, a library of typical chro¬ 
mophores could be used to initialize the class parameters, in¬ 
stead of the proposed class means initialization method. The 
classification process could then perform the function of cor¬ 
recting for uncertainty, errors or local variations in the real op¬ 
tical properties with respect to the prior information. Alterna¬ 
tively, it could be used to label the chromophores found by the 
segmentation process, and identify these as certain tissues such 
as for example ‘oxygenated blood’ or ‘fat’, on the basis of the 
closeness of the recovered means to the characteristic proper¬ 
ties. 

5.7. Additional priors 

In this study we assumed independence between nodal val¬ 
ues, however the mixture of Gaussian model could be used in 
conjunction with a spatial prior. Knowledge of smoothness or 
sparsity properties of the solution could be employed to in¬ 
troduce a homogeneous spatial regularizer such as first-order 
Tikhonov [13] or Total Variation [6, 14]. Knowledge of struc¬ 
tural information, such as that provided by an alternative imag¬ 
ing method or anatomical library, could be exploited by intro¬ 
ducing a spatially varying probability map for the optical prop¬ 
erties. 
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6. Conclusions 


In this paper, we proposed a novel method for perform¬ 
ing image reconstruction in QPAT. We introduced a paramet¬ 
ric class model for the optical parameters, and implemented 
a minimization-based reconstruction algorithm. We suggested 
an automated method by which to initialize the parameters of 
the class model, and proposed a classification algorithm by 
which to progressively update and improve those parameters af¬ 
ter each reconstruction step. We demonstrated though 2D and 
3D numerical examples that the reconstruction-classification 
method allows for the simultaneous recovery of optical absorp¬ 
tion and scattering. In particular, we found that this approach 
delivered superior accuracy in the recovery of the scattering 
with respect to traditional gradient-based reconstruction. 
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